home *** CD-ROM | disk | FTP | other *** search
/ Amiga Collections: Topik / Topik - Disk 39 - Educational (19xx)(Topik Public Domain)(PD)[WB].zip / Topik - Disk 39 - Educational (19xx)(Topik Public Domain)(PD)[WB].adf / Planets / planet_pos.c < prev    next >
C/C++ Source or Header  |  1991-03-09  |  23KB  |  807 lines
  1. /* compute the position of the m`planets and return the epli value */
  2. #include <stdio.h>
  3. #include <math.h>
  4.  
  5. /* crude magnitudes of planets (x100) for output filtering by brightness */
  6.  
  7. #define MAGSOL -9.9
  8. #define MAGMER -1.50
  9. #define MAGVEN -4.0
  10. #define MAGMAR -1.00
  11. #define MAGJUP -2.00
  12. #define MAGSAT 1.00
  13. #define MAGURA 5.90
  14. #define MAGNEP 8.00
  15.  
  16. extern double    rad;
  17. double    htod(), atof(), kepler(), truean();
  18. double    longi(), lati(), poly(), aint(), range();
  19. int    cls(), helio_trans(), geo_trans(), speak();
  20.  
  21. double
  22. planet_pos(jd, T0)
  23. double    jd, T0;
  24. {
  25.     double    l, a, e, i, w1, w2, M, M0, M1, M2, M4, M5, M6, M7, M8;
  26.     double    a0, a1, a2, a3, RA, DEC;
  27.     double    ECC, nu, r, u, ll, b, lonpert, radpert;
  28.     double    esun, Lsun, Cen, Stheta, Snu, Sr;
  29.     double    N, D, epli, thapp, omeg;
  30.     double    nu2, P, Q, S, V, W, ze, l1pert, epert, w1pert, apert;
  31.     double    psi, H, G, eta, th;
  32.  
  33.  
  34.     /* Calculate orbital elements for Mercury */
  35.     a0 = 178.179078;
  36.     a1 = 149474.07078;
  37.     a2 = 0.0003011;
  38.     a3 = 0.0;
  39.     l = poly(a0, a1, a2, a3, T0);
  40.     l = range(l);
  41.     a = 0.3870986;
  42.     a0 = 0.20561421;
  43.     a1 = 0.00002046;
  44.     a2 = -0.000000030;
  45.     e = poly(a0, a1, a2, a3, T0);
  46.     a0 = 7.002881;
  47.     a1 = 0.0018608;
  48.     a2 = -0.0000183;
  49.     i = poly(a0, a1, a2, a3, T0);
  50.     a0 = 28.753753;
  51.     a1 = 0.3702806;
  52.     a2 = 0.0001208;
  53.     w1 = poly(a0, a1, a2, a3, T0);
  54.     w1 = range(w1);
  55.     a0 = 47.145944;
  56.     a1 = 1.1852083;
  57.     a2 = 0.0001739;
  58.     w2 = poly(a0, a1, a2, a3, T0);
  59.     w2 = range(w2);
  60.     M1 = 102.27938 + 149472.51529 * T0 + 0.000007 * T0 * T0;
  61.     M1 = range(M1);
  62.  
  63.     /* solving Kepler find the eccentric anomaly ECC    */
  64.  
  65.     ECC = kepler(e, M1);
  66.     nu = truean(e, ECC);
  67.     r = a * (1.0 - e * cos(ECC * rad));
  68.     u = l + nu - M1 - w2;
  69.     ll = longi(w2, i, u);
  70.     b = lati(u, i);
  71.  
  72.     /* Now make corrections due to perturbations */
  73.     M2 = 212.60322 + 58517.80387 * T0 + 0.001286 * T0 * T0;
  74.     M2 = range(M2);
  75.     M4 = 319.51913 + 19139.85475 * T0 + 0.000181 * T0 * T0;
  76.     M4 = range(M4);
  77.     M5 = 225.32833 + 3034.69202 * T0 - 0.000722 * T0 * T0;
  78.     M5 = range(M5);
  79.     M6 = 175.46622 + 1221.55147 * T0 - 0.000502 * T0 * T0;
  80.     M6 = range(M6);
  81.     lonpert = 0.00204 * cos((5 * M2 - 2 * M1 + 12.220) * rad)
  82.      + 0.00103 * cos((2 * M2 - M1 - 160.692) * rad)
  83.      + 0.00091 * cos((2 * M5 - M1 - 37.003) * rad)
  84.      + 0.00078 * cos((5 * M2 - 3 * M1 + 10.137) * rad);
  85.  
  86.     radpert = 0.000007525 * cos((2 * M5 - M1 + 53.013) * rad)
  87.      + 0.000006802 * cos((5 * M2 - 3 * M1 - 259.918) * rad)
  88.      + 0.000005457 * cos((2 * M2 - 2 * M1 - 71.188) * rad)
  89.      + 0.000003569 * cos((5 * M2 - M1 - 77.75) * rad);
  90.  
  91.     r = r + radpert;
  92.     ll = ll + lonpert;
  93.  
  94.     /* Now find the Longi and radius vector of the sun */
  95.  
  96.     M = 358.47583 + 35999.04975 * T0 - 0.000150 * T0 * T0
  97.         -0.0000033 * T0 * T0 * T0;
  98.     M = range(M);
  99.     esun = 0.01675104 - 0.0000418 * T0 - 0.000000126 * T0 * T0;
  100.     Lsun = 279.69668 + 36000.76892 * T0 + 0.0003025 * T0 * T0;
  101.     Lsun = range(Lsun);
  102.     Cen = (1.919460 - 0.004789 * T0 - 0.000014 * T0 * T0) * sin(M * rad)
  103.      + (0.020094 - 0.000100 * T0) * sin(2 * M * rad)
  104.      + 0.000293 * sin(3 * M * rad);
  105.  
  106.     Stheta = Lsun + Cen;
  107.     Stheta = range(Stheta);
  108.     Snu = M + Cen;
  109.     Sr = 1.0000002 * (1.0 - esun * esun) / (1.0 + esun * cos(Snu * rad));
  110.  
  111.     omeg = 259.18 - 1934.142 * T0;
  112.     thapp = Stheta - 0.00569 - 0.00479 * sin(omeg * rad);
  113.     epli = 23.452294 - 0.0130125 * T0
  114.         -0.00000164 * T0 * T0 + 0.000000503 * T0 * T0 * T0;
  115.  
  116.     N = cos(epli * rad) * sin(thapp * rad);
  117.     D = cos(thapp * rad);
  118.     RA = atan2(N, D) / rad;
  119.     DEC = asin(sin(epli * rad) * sin(thapp * rad)) / rad;
  120.     speak(RA, DEC, Sr, MAGSOL, "PS", "Sol");
  121.     /* tansformation of coordinates on Mercury and output */
  122.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGMER, "PM", "Mercury");
  123.  
  124.     /* Now start on Venus */
  125.     a0 = 342.767053;
  126.     a1 = 58519.21191;
  127.     a2 = 0.0003097;
  128.     a3 = 0.0;
  129.     l = poly(a0, a1, a2, a3, T0);
  130.     l = range(l);
  131.     a = 0.7233316;
  132.     a0 = 0.00682069;
  133.     a1 = -0.00004774;
  134.     a2 = 0.000000091;
  135.     e = poly(a0, a1, a2, a3, T0);
  136.     a0 = 3.393631;
  137.     a1 = 0.0010058;
  138.     a2 = -0.0000010;
  139.     i = poly(a0, a1, a2, a3, T0);
  140.     a0 = 54.384186;
  141.     a1 = 0.5081861;
  142.     a2 = -0.0013864;
  143.     w1 = poly(a0, a1, a2, a3, T0);
  144.     w1 = range(w1);
  145.     a0 = 75.779647;
  146.     a1 = 0.8998500;
  147.     a2 = 0.0004100;
  148.     w2 = poly(a0, a1, a2, a3, T0);
  149.     w2 = range(w2);
  150.     /* Venus has a long period pert that needs to be added before Kelper */
  151.     lonpert = 0.00077 * sin((237.24 + 150.27 * T0) * rad);
  152.     l = l + lonpert;
  153.     M0 = M2 + lonpert;
  154.     ECC = kepler(e, M0);
  155.     nu = truean(e, ECC);
  156.     r = a * (1.0 - e * cos(ECC * rad));
  157.     u = l + nu - M0 - w2;
  158.     ll = longi(w2, i, u);
  159.     b = lati(u, i);
  160.  
  161.     /* now Venus perturbations */
  162.  
  163.     lonpert = 0.00313 * cos((2 * M - 2 * M2 - 148.225) * rad)
  164.      + 0.00198 * cos((3 * M - 3 * M2 + 2.565) * rad)
  165.      + 0.00136 * cos((M - M2 - 119.107) * rad)
  166.      + 0.00096 * cos((3 * M - 2 * M2 - 135.912) * rad)
  167.      + 0.00082 * cos((M5 - M2 - 208.087) * rad);
  168.     ll = ll + lonpert;
  169.  
  170.     radpert = 0.000022501 * cos((2 * M - 2 * M2 - 58.208) * rad)
  171.      + 0.000019045 * cos((3 * M - 3 * M2 + 92.577) * rad)
  172.      + 0.000006887 * cos((M5 - M2 - 118.090) * rad)
  173.      + 0.000005172 * cos((M - M2 - 29.110) * rad)
  174.      + 0.000003620 * cos((5 * M - 4 * M2 - 104.208) * rad)
  175.      + 0.000003283 * cos((4 * M - 4 * M2 + 63.513) * rad)
  176.      + 0.000003074 * cos((2 * M5 - 2 * M2 - 55.167) * rad);
  177.     r = r + radpert;
  178.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGVEN, "PV", "Venus");
  179.  
  180.     /* Now  start the planet Mars */
  181.     a0 = 293.737334;
  182.     a1 = 19141.69551;
  183.     a2 = 0.0003107;
  184.     a3 = 0.0;
  185.     l = poly(a0, a1, a2, a3, T0);
  186.     l = range(l);
  187.     a = 1.5236883;
  188.     a0 = 0.09331290;
  189.     a1 = 0.000092064;
  190.     a2 = -0.000000077;
  191.     e = poly(a0, a1, a2, a3, T0);
  192.     a0 = 1.850333;
  193.     a1 = -0.0006750;
  194.     a2 = 0.0000126;
  195.     i = poly(a0, a1, a2, a3, T0);
  196.     a0 = 285.431761;
  197.     a1 = 1.0697667;
  198.     a2 = 0.0001313;
  199.     a3 = 0.00000414;
  200.     w1 = poly(a0, a1, a2, a3, T0);
  201.     w1 = range(w1);
  202.     a0 = 48.786442;
  203.     a1 = 0.7709917;
  204.     a2 = -0.0000014;
  205.     a3 = -0.00000533;
  206.     w2 = poly(a0, a1, a2, a3, T0);
  207.     w2 = range(w2);
  208.  
  209.     /* Mars has a long period perturbation */
  210.     lonpert = -0.01133 * sin((3 * M5 - 8 * M4 + 4 * M) * rad)
  211.     -0.00933 * cos((3 * M5 - 8 * M4 + 4 * M) * rad);
  212.     l = l + lonpert;
  213.     M0 = M4 + lonpert;
  214.     ECC = kepler(e, M0);
  215.     nu = truean(e, ECC);
  216.     r = a * (1.0 - e * cos(ECC * rad));
  217.     u = l + nu - M0 - w2;
  218.     ll = longi(w2, i, u);
  219.     b = lati(u, i);
  220.  
  221.     /* Now Mars Perturbations */
  222.     lonpert = 0.00705 * cos((M5 - M4 - 48.958) * rad)
  223.      + 0.00607 * cos((2 * M5 - M4 - 188.350) * rad)
  224.      + 0.00445 * cos((2 * M5 - 2 * M4 - 191.897) * rad)
  225.      + 0.00388 * cos((M - 2 * M4 + 20.495) * rad)
  226.      + 0.00238 * cos((M - M4 + 35.097) * rad)
  227.      + 0.00204 * cos((2 * M - 3 * M4 + 158.638) * rad)
  228.      + 0.00177 * cos((3 * M4 - M2 - 57.602) * rad)
  229.      + 0.00136 * cos((2 * M - 4 * M4 + 154.093) * rad)
  230.      + 0.00104 * cos((M5 + 17.618) * rad);
  231.     ll = ll + lonpert;
  232.  
  233.     radpert = 0.000053227 * cos((M5 - M4 + 41.1306) * rad)
  234.      + 0.000050989 * cos((2 * M5 - 2 * M4 - 101.9847) * rad)
  235.      + 0.000038278 * cos((2 * M5 - M4 - 98.3292) * rad)
  236.      + 0.000015996 * cos((M - M4 - 55.555) * rad)
  237.      + 0.000014764 * cos((2 * M - 3 * M4 + 68.622) * rad)
  238.      + 0.000008966 * cos((M5 - 2 * M4 + 43.615) * rad);
  239.     /*
  240.      * broken into two pieces for wimpy C compilers
  241.      */
  242.     radpert += 0.000007914 * cos((3 * M5 - 2 * M4 - 139.737) * rad)
  243.      + 0.000007004 * cos((2 * M5 - 3 * M4 - 102.888) * rad)
  244.      + 0.000006620 * cos((M - 2 * M4 + 113.202) * rad)
  245.      + 0.000004930 * cos((3 * M5 - 3 * M4 - 76.243) * rad)
  246.      + 0.000004693 * cos((3 * M - 5 * M4 + 190.603) * rad)
  247.      + 0.000004571 * cos((2 * M - 4 * M4 + 244.702) * rad)
  248.      + 0.000004409 * cos((3 * M5 - M4 - 115.828) * rad);
  249.     r = r + radpert;
  250.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGMAR, "Pm", "Mars");
  251.  
  252.     /* Now start Jupiter */
  253.     a0 = 238.049257;
  254.     a1 = 3036.301986;
  255.     a2 = 0.0003347;
  256.     a3 = -0.00000165;
  257.     l = poly(a0, a1, a2, a3, T0);
  258.     l = range(l);
  259.     a = 5.202561;
  260.     a0 = 0.04833475;
  261.     a1 = 0.000164180;
  262.     a2 = -0.0000004676;
  263.     a3 = -0.0000000017;
  264.     e = poly(a0, a1, a2, a3, T0);
  265.     a0 = 1.308736;
  266.     a1 = -0.0056961;
  267.     a2 = 0.0000039;
  268.     a3 = 0.0;
  269.     i = poly(a0, a1, a2, a3, T0);
  270.     a0 = 273.277558;
  271.     a1 = 0.5994317;
  272.     a2 = 0.00070405;
  273.     a3 = 0.00000508;
  274.     w1 = poly(a0, a1, a2, a3, T0);
  275.     w1 = range(w1);
  276.     a0 = 99.443414;
  277.     a1 = 1.0105300;
  278.     a2 = 0.00035222;
  279.     a3 = -0.00000851;
  280.     w2 = poly(a0, a1, a2, a3, T0);
  281.     w2 = range(w2);
  282.     /*
  283.      * Now start perturbation calclations
  284.      */
  285.  
  286.     nu2 = T0 / 5.0 + 0.1;
  287.     P = 237.47555 + 3034.9061 * T0;
  288.     Q = 265.91650 + 1222.1139 * T0;
  289.     S = 243.51721 + 428.4677 * T0;
  290.     V = 5.0 * Q - 2.0 * P;
  291.     W = 2.0 * P - 6.0 * Q + 3.0 * S;
  292.     ze = Q - P;
  293.     psi = S - Q;
  294.  
  295.     P = range(P) * rad;
  296.     Q = range(Q) * rad;
  297.     S = range(S) * rad;
  298.     V = range(V) * rad;
  299.     W = range(W) * rad;
  300.     ze = range(ze) * rad;
  301.     psi = range(psi) * rad;
  302.  
  303.     l1pert = (0.331364 - 0.010281 * nu2 - 0.004692 * nu2 * nu2) * sin(V)
  304.      + (0.003228 - 0.064436 * nu2 + 0.002075 * nu2 * nu2) * cos(V)
  305.     -(0.003083 + 0.000275 * nu2 - 0.000489 * nu2 * nu2) * sin(2 * V)
  306.      + 0.002472 * sin(W)
  307.      + 0.013619 * sin(ze)
  308.      + 0.018472 * sin(2 * ze)
  309.      + 0.006717 * sin(3 * ze)
  310.      + 0.002775 * sin(4 * ze)
  311.      + (0.007275 - 0.001253 * nu2) * sin(ze) * sin(Q)
  312.      + 0.006417 * sin(2 * ze) * sin(Q)
  313.      + 0.002439 * sin(3 * ze) * sin(Q);
  314.  
  315.     l1pert = l1pert - (0.033839 + 0.001125 * nu2) * cos(ze) * sin(Q)
  316.     -0.003767 * cos(2 * ze) * sin(Q)
  317.     -(0.035681 + 0.001208 * nu2) * sin(ze) * cos(Q)
  318.     -0.004261 * sin(2 * ze) * cos(Q)
  319.      + 0.002178 * cos(Q)
  320.      + (-0.006333 + 0.001161 * nu2) * cos(ze) * cos(Q)
  321.     -0.006675 * cos(2 * ze) * cos(Q)
  322.     -0.002664 * cos(3 * ze) * cos(Q)
  323.     -0.002572 * sin(ze) * sin(2 * Q)
  324.     -0.003567 * sin(2 * ze) * sin(2 * Q)
  325.      + 0.002094 * cos(ze) * cos(2 * Q)
  326.      + 0.003342 * cos(2 * ze) * cos(2 * Q);
  327.  
  328.     epert = (.0003606 + .0000130 * nu2 - .0000043 * nu2 * nu2) * sin(V)
  329.      + (.0001289 - .0000580 * nu2) * cos(V)
  330.     -.0006764 * sin(ze) * sin(Q)
  331.     -.0001110 * sin(2 * ze) * sin(Q)
  332.     -.0000224 * sin(3 * ze) * sin(Q)
  333.     -.0000204 * sin(Q)
  334.      + (.0001284 + .0000116 * nu2) * cos(ze) * sin(Q)
  335.      + .0000188 * cos(2 * ze) * sin(Q)
  336.      + (.0001460 + .0000130 * nu2) * sin(ze) * cos(Q)
  337.      + .0000224 * sin(2 * ze) * cos(Q)
  338.     -.0000817 * cos(Q);
  339.  
  340.     epert = epert + .0006074 * cos(ze) * cos(Q)
  341.      + .0000992 * cos(2 * ze) * cos(Q)
  342.      + .0000508 * cos(3 * ze) * cos(Q)
  343.      + .0000230 * cos(4 * ze) * cos(Q)
  344.      + .0000108 * cos(5 * ze) * cos(Q)
  345.     -(.0000956 + .0000073 * nu2) * sin(ze) * sin(2 * Q)
  346.      + .0000448 * sin(2 * ze) * sin(2 * Q)
  347.      + .0000137 * sin(3 * ze) * sin(2 * Q)
  348.      + (-.0000997 + .0000108 * nu2) * cos(ze) * sin(2 * Q)
  349.      + .0000480 * cos(2 * ze) * sin(2 * Q);
  350.  
  351.     epert = epert + .0000148 * cos(3 * ze) * sin(2 * Q)
  352.      + (-.0000956 + .0000099 * nu2) * sin(ze) * cos(2 * Q)
  353.      + .0000490 * sin(2 * ze) * cos(2 * Q)
  354.      + .0000158 * sin(3 * ze) * cos(2 * Q)
  355.      + .0000179 * cos(2 * Q)
  356.      + (.0001024 + .0000075 * nu2) * cos(ze) * cos(2 * Q)
  357.     -.0000437 * cos(2 * ze) * cos(2 * Q)
  358.     -.0000132 * cos(3 * ze) * cos(2 * Q);
  359.  
  360.     w1pert = (0.007192 - 0.003147 * nu2) * sin(V)
  361.      + (-0.020428 - 0.000675 * nu2 + 0.000197 * nu2 * nu2) * cos(V)
  362.      + (0.007269 + 0.000672 * nu2) * sin(ze) * sin(Q)
  363.     -0.004344 * sin(Q)
  364.      + 0.034036 * cos(ze) * sin(Q)
  365.      + 0.005614 * cos(2 * ze) * sin(Q)
  366.      + 0.002964 * cos(3 * ze) * sin(Q)
  367.      + 0.037761 * sin(ze) * cos(Q);
  368.  
  369.     w1pert = w1pert
  370.          + 0.006158 * sin(2 * ze) * cos(Q)
  371.     -0.006603 * cos(ze) * cos(Q)
  372.     -0.005356 * sin(ze) * sin(2 * Q)
  373.      + 0.002722 * sin(2 * ze) * sin(2 * Q)
  374.      + 0.004483 * cos(ze) * sin(2 * Q)
  375.     -0.002642 * cos(2 * ze) * sin(2 * Q)
  376.      + 0.004403 * sin(ze) * cos(2 * Q)
  377.     -0.002536 * sin(2 * ze) * cos(2 * Q)
  378.      + 0.005547 * cos(ze) * cos(2 * Q)
  379.     -0.002689 * cos(2 * ze) * cos(2 * Q);
  380.  
  381.     lonpert = l1pert - (w1pert / e);
  382.     l = l + l1pert;
  383.     M0 = M5 + lonpert;
  384.     e = e + epert;
  385.     w1 = w1 + w1pert;
  386.  
  387.     apert = -.000263 * cos(V)
  388.      + .000205 * cos(ze)
  389.      + .000693 * cos(2 * ze)
  390.      + .000312 * cos(3 * ze)
  391.      + .000147 * cos(4 * ze)
  392.      + .000299 * sin(ze) * sin(Q)
  393.      + .000181 * cos(2 * ze) * sin(Q)
  394.      + .000204 * sin(2 * ze) * cos(Q)
  395.      + .000111 * sin(3 * ze) * cos(Q)
  396.     -.000337 * cos(ze) * cos(Q)
  397.     -.000111 * cos(2 * ze) * cos(Q);
  398.  
  399.     a = a + apert;
  400.     ECC = kepler(e, M0);
  401.     nu = truean(e, ECC);
  402.     r = a * (1.0 - e * cos(ECC * rad));
  403.     u = l + nu - M0 - w2;
  404.     ll = longi(w2, i, u);
  405.     b = lati(u, i);
  406.  
  407.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGJUP, "PJ", "Jupiter");
  408.  
  409.     /* Now start Saturn */
  410.  
  411.     a0 = 266.564377;
  412.     a1 = 1223.509884;
  413.     a2 = 0.0003245;
  414.     a3 = -0.0000058;
  415.     l = poly(a0, a1, a2, a3, T0);
  416.     l = range(l);
  417.     a = 9.554747;
  418.     a0 = 0.05589232;
  419.     a1 = -0.00034550;
  420.     a2 = -0.000000728;
  421.     a3 = 0.00000000074;
  422.     e = poly(a0, a1, a2, a3, T0);
  423.     a0 = 2.492519;
  424.     a1 = -0.0039189;
  425.     a2 = -0.00001549;
  426.     a3 = 0.00000004;
  427.     i = poly(a0, a1, a2, a3, T0);
  428.     a0 = 338.307800;
  429.     a1 = 1.0852207;
  430.     a2 = 0.00097854;
  431.     a3 = 0.00000992;
  432.     w1 = poly(a0, a1, a2, a3, T0);
  433.     w1 = range(w1);
  434.     a0 = 112.790414;
  435.     a1 = 0.8731951;
  436.     a2 = -0.00015218;
  437.     a3 = -0.00000531;
  438.     w2 = poly(a0, a1, a2, a3, T0);
  439.     w2 = range(w2);
  440.  
  441.     /* Now Saturn's perturbations */
  442.  
  443.     l1pert = (-0.814181 + 0.018150 * nu2 + 0.016714 * nu2 * nu2) * sin(V)
  444.      + (-0.010497 + 0.160906 * nu2 - 0.004100 * nu2 * nu2) * cos(V)
  445.      + 0.007581 * sin(2 * V)
  446.     -0.007986 * sin(W)
  447.     -0.148811 * sin(ze)
  448.     -0.040786 * sin(2 * ze)
  449.     -0.015208 * sin(3 * ze)
  450.     -0.006339 * sin(4 * ze)
  451.     -0.006244 * sin(Q);
  452.     l1pert = l1pert
  453.          + (0.008931 + 0.002728 * nu2) * sin(ze) * sin(Q)
  454.     -0.016500 * sin(2 * ze) * sin(Q)
  455.     -0.005775 * sin(3 * ze) * sin(Q)
  456.      + (0.081344 + 0.003206 * nu2) * cos(ze) * sin(Q)
  457.      + 0.015019 * cos(2 * ze) * sin(Q)
  458.      + (0.085581 + 0.002494 * nu2) * sin(ze) * cos(Q)
  459.      + (0.025328 - 0.003117 * nu2) * cos(ze) * cos(Q);
  460.     l1pert = l1pert
  461.          + 0.014394 * cos(2 * ze) * cos(Q)
  462.      + 0.006319 * cos(3 * ze) * cos(Q)
  463.      + 0.006369 * sin(ze) * sin(2 * Q)
  464.      + 0.009156 * sin(2 * ze) * sin(2 * Q)
  465.      + 0.007525 * sin(3 * psi) * sin(2 * Q)
  466.     -0.005236 * cos(ze) * cos(2 * Q)
  467.     -0.007736 * cos(2 * ze) * cos(2 * Q)
  468.     -0.007528 * cos(3 * psi) * cos(2 * Q);
  469.  
  470.     epert = (-.0007927 + .0002548 * nu2 + .0000091 * nu2 * nu2) * sin(V)
  471.      + (.0013381 + .0001226 * nu2 - .0000253 * nu2 * nu2) * cos(V)
  472.      + (.0000248 - .0000121 * nu2) * sin(2 * V)
  473.     -(.0000305 + .0000091 * nu2) * cos(2 * V)
  474.      + .0000412 * sin(2 * ze)
  475.      + .0012415 * sin(Q)
  476.      + (.0000390 - .0000617 * nu2) * sin(ze) * sin(Q)
  477.      + (.0000165 - .0000204 * nu2) * sin(2 * ze) * sin(Q)
  478.      + .0026599 * cos(ze) * sin(Q)
  479.     -.0004687 * cos(2 * ze) * sin(Q);
  480.     epert = epert
  481.         -.0001870 * cos(3 * ze) * sin(Q)
  482.     -.0000821 * cos(4 * ze) * sin(Q)
  483.     -.0000377 * cos(5 * ze) * sin(Q)
  484.      + .0000497 * cos(2 * psi) * sin(Q)
  485.      + (.0000163 - .0000611 * nu2) * cos(Q)
  486.     -.0012696 * sin(ze) * cos(Q)
  487.     -.0004200 * sin(2 * ze) * cos(Q)
  488.     -.0001503 * sin(3 * ze) * cos(Q)
  489.     -.0000619 * sin(4 * ze) * cos(Q)
  490.     -.0000268 * sin(5 * ze) * cos(Q);
  491.     epert = epert
  492.         -(.0000282 + .0001306 * nu2) * cos(ze) * cos(Q)
  493.      + (-.0000086 + .0000230 * nu2) * cos(2 * ze) * cos(Q)
  494.      + .0000461 * sin(2 * psi) * cos(Q)
  495.     -.0000350 * sin(2 * Q)
  496.      + (.0002211 - .0000286 * nu2) * sin(ze) * sin(2 * Q)
  497.     -.0002208 * sin(2 * ze) * sin(2 * Q)
  498.     -.0000568 * sin(3 * ze) * sin(2 * Q)
  499.     -.0000346 * sin(4 * ze) * sin(2 * Q)
  500.     -(.0002780 + .0000222 * nu2) * cos(ze) * sin(2 * Q)
  501.      + (.0002022 + .0000263 * nu2) * cos(2 * ze) * sin(2 * Q);
  502.     epert = epert
  503.          + .0000248 * cos(3 * ze) * sin(2 * Q)
  504.      + .0000242 * sin(3 * psi) * sin(2 * Q)
  505.      + .0000467 * cos(3 * psi) * sin(2 * Q)
  506.     -.0000490 * cos(2 * Q)
  507.     -(.0002842 + .0000279 * nu2) * sin(ze) * cos(2 * Q)
  508.      + (.0000128 + .0000226 * nu2) * sin(2 * ze) * cos(2 * Q)
  509.      + .0000224 * sin(3 * ze) * cos(2 * Q)
  510.      + (-.0001594 + .0000282 * nu2) * cos(ze) * cos(2 * Q)
  511.      + (.0002162 - .0000207 * nu2) * cos(2 * ze) * cos(2 * Q)
  512.      + .0000561 * cos(3 * ze) * cos(2 * Q);
  513.     epert = epert
  514.          + .0000343 * cos(4 * ze) * cos(2 * Q)
  515.      + .0000469 * sin(3 * psi) * cos(2 * Q)
  516.     -.0000242 * cos(3 * psi) * cos(2 * Q)
  517.     -.0000205 * sin(ze) * sin(3 * Q)
  518.      + .0000262 * sin(3 * ze) * sin(3 * Q)
  519.      + .0000208 * cos(ze) * cos(3 * Q)
  520.     -.0000271 * cos(3 * ze) * cos(3 * Q)
  521.     -.0000382 * cos(3 * ze) * sin(4 * Q)
  522.     -.0000376 * sin(3 * ze) * cos(4 * Q);
  523.  
  524.     w1pert = (0.077108 + 0.007186 * nu2 - 0.001533 * nu2 * nu2) * sin(V)
  525.      + (0.045803 - 0.014766 * nu2 - 0.000536 * nu2 * nu2) * cos(V)
  526.     -0.007075 * sin(ze)
  527.     -0.075825 * sin(ze) * sin(Q)
  528.     -0.024839 * sin(2 * ze) * sin(Q)
  529.     -0.008631 * sin(3 * ze) * sin(Q)
  530.     -0.072586 * cos(Q)
  531.     -0.150383 * cos(ze) * cos(Q)
  532.      + 0.026897 * cos(2 * ze) * cos(Q)
  533.      + 0.010053 * cos(3 * ze) * cos(Q);
  534.     w1pert = w1pert
  535.         -(0.013597 + 0.001719 * nu2) * sin(ze) * sin(2 * Q)
  536.      + (-0.007742 + 0.001517 * nu2) * cos(ze) * sin(2 * Q)
  537.      + (0.013586 - 0.001375 * nu2) * cos(2 * ze) * sin(2 * Q)
  538.      + (-0.013667 + 0.001239 * nu2) * sin(ze) * cos(2 * Q)
  539.      + 0.011981 * sin(2 * ze) * cos(2 * Q)
  540.      + (0.014861 + 0.001136 * nu2) * cos(ze) * cos(2 * Q)
  541.     -(0.013064 + 0.001628 * nu2) * cos(2 * ze) * cos(2 * Q);
  542.  
  543.     lonpert = l1pert - (w1pert / e);
  544.     l = l + l1pert;
  545.     M0 = M6 + lonpert;
  546.     e = e + epert;
  547.     w1 = w1 + w1pert;
  548.  
  549.     apert = .000572 * sin(V) - .001590 * sin(2 * ze) * cos(Q)
  550.      + .002933 * cos(V) - .000647 * sin(3 * ze) * cos(Q)
  551.      + .033629 * cos(ze) - .000344 * sin(4 * ze) * cos(Q)
  552.     -.003081 * cos(2 * ze) + .002885 * cos(ze) * cos(Q)
  553.     -.001423 * cos(3 * ze) + (.002172 + .000102 * nu2) * cos(2 * ze) * cos(Q)
  554.     -.000671 * cos(4 * ze) + .000296 * cos(3 * ze) * cos(Q)
  555.     -.000320 * cos(5 * ze) - .000267 * sin(2 * ze) * sin(2 * Q);
  556.     apert = apert
  557.          + .001098 * sin(Q) - .000778 * cos(ze) * sin(2 * Q)
  558.     -.002812 * sin(ze) * sin(Q) + .000495 * cos(2 * ze) * sin(2 * Q)
  559.      + .000688 * sin(2 * ze) * sin(Q) + .000250 * cos(3 * ze) * sin(2 * Q);
  560.     apert = apert
  561.         -.000393 * sin(3 * ze) * sin(Q)
  562.     -.000228 * sin(4 * ze) * sin(Q)
  563.      + .002138 * cos(ze) * sin(Q)
  564.     -.000999 * cos(2 * ze) * sin(Q)
  565.     -.000642 * cos(3 * ze) * sin(Q)
  566.     -.000325 * cos(4 * ze) * sin(Q)
  567.     -.000890 * cos(Q)
  568.      + .002206 * sin(ze) * cos(Q);
  569.     apert = apert
  570.         -.000856 * sin(ze) * cos(2 * Q)
  571.      + .000441 * sin(2 * ze) * cos(2 * Q)
  572.      + .000296 * cos(2 * ze) * cos(2 * Q)
  573.      + .000211 * cos(3 * ze) * cos(2 * Q)
  574.     -.000427 * sin(ze) * sin(3 * Q)
  575.      + .000398 * sin(3 * ze) * sin(3 * Q)
  576.      + .000344 * cos(ze) * cos(3 * Q)
  577.     -.000427 * cos(3 * ze) * cos(3 * Q);
  578.  
  579.     a = a + apert;
  580.     ECC = kepler(e, M0);
  581.     nu = truean(e, ECC);
  582.     r = a * (1.0 - e * cos(ECC * rad));
  583.     u = l + nu - M0 - w2;
  584.     ll = longi(w2, i, u);
  585.     b = lati(u, i);
  586.  
  587.     b = b
  588.          + 0.000747 * cos(ze) * sin(Q)
  589.      + 0.001069 * cos(ze) * cos(Q)
  590.      + 0.002108 * sin(2 * ze) * sin(2 * Q)
  591.      + 0.001261 * cos(2 * ze) * sin(2 * Q)
  592.      + 0.001236 * sin(2 * ze) * cos(2 * Q)
  593.     -0.002075 * cos(2 * ze) * cos(2 * Q);
  594.  
  595.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGSAT, "Ps", "Saturn");
  596.  
  597.     /* Now Start on Uranus */
  598.     a0 = 244.197470;
  599.     a1 = 429.863546;
  600.     a2 = 0.0003160;
  601.     a3 = -0.00000060;
  602.     l = poly(a0, a1, a2, a3, T0);
  603.     l = range(l);
  604.     a = 19.21814;
  605.     a0 = 0.0463444;
  606.     a1 = -0.00002658;
  607.     a2 = 0.000000077;
  608.     a3 = 0.0;
  609.     e = poly(a0, a1, a2, a3, T0);
  610.     a0 = 0.772464;
  611.     a1 = 0.0006253;
  612.     a2 = 0.0000395;
  613.     i = poly(a0, a1, a2, a3, T0);
  614.     a0 = 98.071581;
  615.     a1 = 0.9857650;
  616.     a2 = -0.0010745;
  617.     a3 = -0.00000061;
  618.     w1 = poly(a0, a1, a2, a3, T0);
  619.     w1 = range(w1);
  620.     a0 = 73.477111;
  621.     a1 = 0.4986678;
  622.     a2 = 0.0013117;
  623.     a3 = 0.0;
  624.     w2 = poly(a0, a1, a2, a3, T0);
  625.     w2 = range(w2);
  626.     M7 = 72.64878 + 428.37911 * T0 + 0.000079 * T0 * T0;
  627.     M7 = range(M7);
  628.     /* now perturbation corrections for Uranus */
  629.     G = 83.76922 + 218.4901 * T0;
  630.     S = S / rad;
  631.     P = P / rad;
  632.     Q = Q / rad;
  633.     H = 2.0 * G - S;
  634.     ze = S - P;
  635.     eta = S - Q;
  636.     th = G - S;
  637.     S = S * rad;
  638.     G = range(G) * rad;
  639.     P = P * rad;
  640.     Q = Q * rad;
  641.     H = range(H) * rad;
  642.     ze = range(ze) * rad;
  643.     eta = range(eta) * rad;
  644.     th = range(th) * rad;
  645.  
  646.     l1pert = (0.864319 - 0.001583 * nu2) * sin(H)
  647.      + (0.082222 - 0.006833 * nu2) * cos(H)
  648.      + 0.036017 * sin(2 * H)
  649.     -0.003019 * cos(2 * H)
  650.      + 0.008122 * sin(W);
  651.  
  652.     epert = (-.0003349 + .0000163 * nu2) * sin(H)
  653.      + .0020981 * cos(H)
  654.      + .0001311 * cos(H);
  655.  
  656.     w1pert = 0.120303 * sin(H)
  657.      + (0.019472 - 0.000947 * nu2) * cos(H)
  658.      + 0.006197 * sin(2 * H);
  659.  
  660.     lonpert = l1pert - (w1pert / e);
  661.     l = l + l1pert;
  662.     M0 = M7 + lonpert;
  663.     e = e + epert;
  664.     w1 = w1 + w1pert;
  665.  
  666.     a = a - 0.003825 * cos(H);
  667.     ECC = kepler(e, M0);
  668.     nu = truean(e, ECC);
  669.     r = a * (1.0 - e * cos(ECC * rad));
  670.     u = l + nu - M0 - w2;
  671.     ll = longi(w2, i, u);
  672.     b = lati(u, i);
  673.  
  674.     ll = ll
  675.          + (0.010122 - 0.000988 * nu2) * sin(S + eta)
  676.      + (-0.038581 + 0.002031 * nu2 - 0.001910 * nu2 * nu2) * cos(S + eta)
  677.      + (0.034964 - 0.001038 * nu2 + 0.000868 * nu2 * nu2) * cos(2 * S + eta)
  678.      + 0.005594 * sin(S + 3 * th);
  679.     ll = ll
  680.         -0.014808 * sin(ze)
  681.     -0.005794 * sin(eta)
  682.      + 0.002347 * cos(eta)
  683.      + 0.009872 * sin(th)
  684.      + 0.008803 * sin(2 * th)
  685.     -0.004308 * sin(3 * th);
  686.  
  687.     b = b
  688.          + (0.000458 * sin(eta) - 0.000642 * cos(eta) - 0.000517 * cos(4 * th))
  689.     *sin(S)
  690.     -(0.000347 * sin(eta) + 0.000853 * cos(eta) + 0.000517 * sin(4 * eta))
  691.     *cos(S)
  692.      + 0.000403 * (cos(2 * th) * sin(2 * S) + sin(2 * th) * cos(2 * S));
  693.  
  694.     r = r
  695.         -.025948
  696.          + .004985 * cos(ze)
  697.     -.001230 * cos(S)
  698.      + .003354 * cos(eta)
  699.      + (.005795 * cos(S) - .001165 * sin(S) + .001388 * cos(2 * S)) * sin(eta)
  700.      + (.001351 * cos(S) + .005702 * sin(S) + .001388 * sin(2 * S)) * cos(eta)
  701.      + .000904 * cos(2 * th)
  702.      + .000894 * (cos(th) - cos(3 * th));
  703.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGURA, "PU", "Uranus");
  704.  
  705.     /* Now start Neptune */
  706.     a0 = 84.457994;
  707.     a1 = 219.885914;
  708.     a2 = 0.0003205;
  709.     a3 = -0.00000060;
  710.     l = poly(a0, a1, a2, a3, T0);
  711.     l = range(l);
  712.     a = 30.10957;
  713.     a0 = 0.00899704;
  714.     a1 = 0.000006330;
  715.     a2 = -0.000000002;
  716.     a3 = 0.0;
  717.     e = poly(a0, a1, a2, a3, T0);
  718.     a0 = 1.779242;
  719.     a1 = -0.0095436;
  720.     a2 = -0.0000091;
  721.     i = poly(a0, a1, a2, a3, T0);
  722.     a0 = 276.045975;
  723.     a1 = 0.3256394;
  724.     a2 = 0.00014095;
  725.     a3 = 0.000004113;
  726.     w1 = poly(a0, a1, a2, a3, T0);
  727.     w1 = range(w1);
  728.     a0 = 130.681389;
  729.     a1 = 1.0989350;
  730.     a2 = 0.00024987;
  731.     a3 = -0.000004718;
  732.     w2 = poly(a0, a1, a2, a3, T0);
  733.     w2 = range(w2);
  734.     M8 = 37.73063 + 218.46134 * T0 - 0.000070 * T0 * T0;
  735.  
  736.     /* now perturbation corrections for neptune */
  737.     G = G / rad;
  738.     P = P / rad;
  739.     Q = Q / rad;
  740.     S = S / rad;
  741.     ze = G - P;
  742.     eta = G - Q;
  743.     th = G - S;
  744.     G = G * rad;
  745.     P = P * rad;
  746.     Q = Q * rad;
  747.     S = S * rad;
  748.     ze = range(ze) * rad;
  749.     eta = range(eta) * rad;
  750.     th = range(th) * rad;
  751.  
  752.     l1pert = (-0.589833 + 0.001089 * nu2) * sin(H)
  753.      + (-0.056094 + 0.004658 * nu2) * cos(H)
  754.     -0.024286 * sin(2 * H);
  755.  
  756.     epert = .0004389 * sin(H)
  757.      + .0004262 * cos(H)
  758.      + .0001129 * sin(2 * H)
  759.      + .0001089 * cos(2 * H);
  760.  
  761.     w1pert = 0.024039 * sin(H)
  762.     -0.025303 * cos(H)
  763.      + 0.006206 * sin(2 * H)
  764.     -0.005992 * cos(2 * H);
  765.  
  766.  
  767.     lonpert = l1pert - (w1pert / e);
  768.     l = l + l1pert;
  769.     M0 = M8 + lonpert;
  770.     e = e + epert;
  771.     w1 = w1 + w1pert;
  772.  
  773.     a = a - 0.000817 * sin(H)
  774.      + 0.008189 * cos(H)
  775.      + 0.000781 * cos(2 * H);
  776.  
  777.     ECC = kepler(e, M0);
  778.     nu = truean(e, ECC);
  779.     r = a * (1.0 - e * cos(ECC * rad));
  780.     u = l + nu - M0 - w2;
  781.     ll = longi(w2, i, u);
  782.     b = lati(u, i);
  783.  
  784.     ll = ll
  785.         -0.009556 * sin(ze)
  786.     -0.005178 * sin(eta)
  787.      + 0.002572 * sin(2 * th)
  788.     -0.002972 * cos(2 * th) * sin(G)
  789.     -0.002833 * sin(2 * th) * cos(G);
  790.  
  791.     b = b
  792.          + 0.000336 * cos(2 * th) * sin(G)
  793.      + 0.000364 * sin(2 * th) * cos(G);
  794.  
  795.     r = r
  796.         -.040596
  797.          + .004992 * cos(ze)
  798.      + .002744 * cos(eta)
  799.      + .002044 * cos(th)
  800.      + .001051 * cos(2 * th);
  801.     helio_trans(r, b, ll, Stheta, Sr, epli, MAGNEP, "PN", "Neptune");
  802.  
  803.     return epli;
  804. }
  805.  
  806.  
  807.